In [ ]:
import eex
import os
import pandas as pd
import numpy as np
In [ ]:
# Create empty data layer
dl = eex.datalayer.DataLayer("butane", backend="Memory")
dl.summary()
In [ ]:
"""
First, we add atoms to the system. Atoms have associated metadata. The possible atom metadata is listed here.
"""
dl.list_valid_atom_properties()
In [ ]:
"""
TOPOLOGY:
Information can be added to the datalayer in the form of pandas dataframes. Here, we add atom metadata.
The name of the column corresponds to the atom property.
Populate empty dataframe with relevant information and add to EEX datalayer
"""
# Create empty dataframe
atom_df = pd.DataFrame()
# Create atomic system using pandas dataframe
atom_df["atom_index"] = np.arange(0,4)
atom_df["molecule_index"] = [int(x) for x in np.zeros(4)]
atom_df["residue_index"] = [int(x) for x in np.zeros(4)]
atom_df["atom_name"] = ["C1", "C2", "C3", "C4"]
atom_df["charge"] = np.zeros(4)
atom_df["atom_type"] = [1, 2, 2, 1]
atom_df["X"] = [0, 0, 0, -1.474]
atom_df["Y"] = [-0.4597, 0, 1.598, 1.573]
atom_df["Z"] = [-1.5302, 0, 0, -0.6167]
atom_df["mass"] = [15.0452, 14.02658, 14.02658, 15.0452]
# Add atoms to datalayer
dl.add_atoms(atom_df, by_value=True)
# Print datalayer information
dl.summary()
# Print stored atom properties
dl.get_atoms(properties=None, by_value=True)
In [ ]:
"""
TOPOLOGY:
The EEX datalayer now contains four nonbonded atoms. To create butane, atoms must be bonded
to one another.
Add bonds to system
"""
# Create empty dataframes for bonds
bond_df = pd.DataFrame()
# Create column names. Here, "term_index" refers to the bond type index.
# i.e. - if all bonds are the same type, they will have the same term index
bond_column_names = ["atom1", "atom2", "term_index"]
# Create corresponding data. The first row specifies that atom0 is bonded
# to atom 1 and has bond_type id 0
bond_data = np.array([[0, 1, 0,],
[1, 2, 0],
[2, 3, 0]])
for num, name in enumerate(bond_column_names):
bond_df[name] = bond_data[:,num]
dl.add_bonds(bond_df)
In [ ]:
dl.summary()
In [ ]:
"""
TOPOLOGY:
Add angles and dihedrals to system.
"""
# Follow similar procedure as for bonds
angle_df = pd.DataFrame()
dihedral_df = pd.DataFrame()
angle_column_names = ["atom1", "atom2", "atom3", "term_index"]
dihedral_column_names = ["atom1", "atom2", "atom3", "atom4", "term_index"]
angle_data = np.array([[0, 1, 2, 0,],
[1, 2, 3, 0],])
dihedral_data = np.array([[0, 1, 2, 3, 0,]])
for num, name in enumerate(angle_column_names):
angle_df[name] = angle_data[:,num]
dl.add_angles(angle_df)
for num, name in enumerate(dihedral_column_names):
dihedral_df[name] = dihedral_data[:,num]
dl.add_dihedrals(dihedral_df)
In [ ]:
dl.summary()
So far, only the topology and coordinates of the system are specified, and we are not able to calculate an energy.
To calculate the energy, we need to define the functional form of bond, angle, dihedral, and nonbonded interactions and the associated constants.
In this demo, we store the parameters for the TraPPE United Atom forcefield with harmonic bonds.
\begin{equation*} \ U_{total} = \sum_{bonds}{k_{b}(r-r_{0})^2} + \sum_{angles}{k_{\theta} (\theta - \theta_{eq} )^{2}} + \sum_{dihedrals}{c_{1}[1 + cos(\phi)] + c_{2}[1 - cos(2\phi)] + c_{3}[1 + cos(3\phi)]} + \sum_{i=1}^{N-1}{\sum_{j=i+1}^{N}{ 4\epsilon_{ij}[(\frac{\sigma_{ij}}r_{ij})^{12} - (\frac{\sigma_{ij}}r_{ij})^6] }} \end{equation*}
In [ ]:
"""
EEX FORCE FIELD PARAMETERS
A main component of EEX is internally stored metadata which defines details functional forms including form, constants,
unit types, and default units (if the user does not overrride this option).
This metadata is stored as human readable dictionaries which can easily be added or maninpulated.
"""
# Here, we examine the metadata present in the bond metadata for a harmonic bond
bond_metadata = eex.metadata.two_body_terms.two_body_metadata
for k, v in bond_metadata["forms"]["harmonic"].items():
print(k, v)
In [ ]:
"""
FORCE FIELD PARAMETERS
To add bonds (or other parameters) using this metadata, the user specifies the form using a keyword ("harmonic") that
matches to EEX's metadata.
Values for the contstants are passed using a dictionary with the 'parameters' defined in the metadata as keys.
Each bond type is given a uid, and default dimensions may be overwritten by the user using a dictionary
and the 'utype' argument
"""
# Here, in add_term_parameter, the first argument is the term order. '2'
# corresponds to bonded atoms.
dl.add_term_parameter(2, "harmonic", {'K': 300.9, 'R0': 1.540}, uid=0, utype={'K':"kcal * mol **-1 * angstrom ** -2",
'R0': "angstrom"})
In [ ]:
# If units or parameters are not compatible with the metadata, the datalayer will not allow storage of the parameter.
# Here, we have changed 'K' to simply "kcal". This will fail (uncomment to test)
#dl.add_term_parameter(2, "harmonic", {'K': 300.9, 'R0': 1.540}, uid=0, utype={'K':"kcal",'R0': "angstrom"})
In [ ]:
## Add harmonic angle parameters
dl.add_term_parameter(3, "harmonic", {'K': 62.100, 'theta0': 114}, uid=0, utype={'K':'kcal * mol ** -1 * radian ** -2',
'theta0': 'degree'})
# Add OPLS dihedral parameter
dl.add_term_parameter(4, "opls", {'K_1': 1.41103414, 'K_2': -0.27101489,
'K_3': 3.14502869, 'K_4': 0}, uid=0, utype={'K_1': 'kcal * mol ** -1',
'K_2': 'kcal * mol ** -1',
'K_3': 'kcal * mol ** -1',
'K_4': 'kcal * mol ** -1'})
In [ ]:
"""
NONBOND PARAMETERS
For nonbond parametets, we currently provide support for Lennard Jones and Buckingham potentials
Most programs use pair-wise Lennard Jones potentials for nonbond interactions. Our internal metadata stores these as A
and B parameters. However, uses may specify other forms such as epsilon/sigma, epsilon, Rmin, etc.
Lennard Jones parameters can be added as a pair (atom_type1, atom_type2) or for a single atom type with a mixing rule.
"""
dl.add_nb_parameter(atom_type=1, nb_name="LJ",
nb_model="epsilon/sigma", nb_parameters={'sigma': 3.75, 'epsilon': 0.1947460018},
utype={'sigma': 'angstrom', 'epsilon': 'kcal * mol ** -1'})
dl.add_nb_parameter(atom_type=2, nb_name="LJ",
nb_model="epsilon/sigma", nb_parameters={'sigma': 3.95, 'epsilon': 0.0914112887},
utype={'sigma': 'angstrom', 'epsilon': 'kcal * mol ** -1'})
dl.set_mixing_rule('lorentz-berthelot')
In [ ]:
# Retrieve stored parameters
print("All stored parameters\n", dl.list_nb_parameters("LJ"), "\n\n")
# To apply the mixing rule:
dl.build_LJ_mixing_table()
print("All stored parameters\n", dl.list_nb_parameters("LJ"), "\n\n")
# These can also be retrieved for only single atoms, or for atom pairs by using itype='single' or itype='pairs'
pair_interactions = dl.list_nb_parameters("LJ", itype="pair")
print("Pair parameters\n", pair_interactions)
Alternatively, these could have been set directly as pairs without a mixing rule.
# Add NB parameters with pairs
dl.add_nb_parameter(atom_type=1, atom_type2=1, nb_name="LJ", nb_model="AB", nb_parameters=[1.0, 1.0])
dl.add_nb_parameter(atom_type=1, atom_type2=2, nb_name="LJ", nb_model="epsilon/sigma", nb_parameters=[1.0, 1.0])
dl.add_nb_parameter(atom_type=2, atom_type2=2, nb_name="LJ", nb_model="epsilon/sigma", nb_parameters=[1.0, 1.0])
In [ ]:
dl.summary()
In [ ]:
# Evaluate system energy
energy_system1 = dl.evaluate(utype="kcal * mol ** -1")
print(energy_system1)
Typically, this storage would not be done by hand as shown above. Instead, readers and writers for specific softwares are used.
Below, we use the plugin for AMBER with EEX to read in an amber file with information for the a butane molecule which is equivalent to the one created in the first datalayer. The EEX translator uses the functions displayed above to store all information from the amber prmtop and inpcrd files in the datalayer.
In [ ]:
# Preview an amber prmtop (parameter-topology file) for Amber.
butane_file = os.path.join("..", "examples", "amber","alkanes", "trappe_butane_single_molecule.prmtop")
f = open(butane_file)
print(f.read())
f.close()
In [ ]:
# Create new datalayer and populate using amber reader
dl_amber = eex.datalayer.DataLayer("butane_amber")
eex.translators.amber.read_amber_file(dl_amber, butane_file)
dl_amber.summary()
The summary shows a system with 4 atom, 3 bonds, 2 angles and 4 dihedrals. This differs from the first datalayer in the number of dihedrals and the number of dihedral parameters. However, the evaluated system energy is equivalent.
This is because AMBER stores dihedral angles using a different functional form. Instead of a single equation, dihedrals with multiple terms are built from multiple harmonic equations. The equations are equivalent when evaluated.
\begin{equation*} \sum_{dihedrals}{(V_{0}[1 + cos(0\phi)]}) + \sum_{dihedrals}{(V_{1}[1 + cos(n\phi)]}) + \sum_{dihedrals}{(V_{2}[1 + cos(2\phi - \pi)]}) + \sum_{dihedrals}{(V_{3}[1 + cos(3\phi )]}) = \sum_{dihedrals}{c_{0} + c_{1}[1 + cos(\phi)] + c_{2}[1 - cos(2\phi)] + c_{3}[1 + cos(3\phi)]} \end{equation*}Although not not implemented yet, EEX should eventually be able to identify and perform a translation between equivalent functional forms.
In [ ]:
energy_system2 = dl_amber.evaluate(utype="kcal * mol ** -1")
In [ ]:
for k in energy_system1:
energy_difference = energy_system1[k] - energy_system2[k]
print(k," difference:\t %.3f" % energy_difference)
In [ ]:
# Compare stored NB parameters
eex.testing.dict_compare(dl_amber.list_nb_parameters("LJ"), dl.list_nb_parameters("LJ", itype="pair"))
In [ ]:
# We can now write the amber file we read for lammps.
eex.translators.lammps.write_lammps_file(dl_amber, "output_lammps.data", "output_lammps.in")
# Write a local copy of the amber datalayer for amber.
eex.translators.amber.write_amber_file(dl_amber, "amber_output.prmtop")
In [ ]:
## Read the written file into a datalayer ##
dl_lammps = eex.datalayer.DataLayer("butane_lammps")
eex.translators.lammps.read_lammps_input_file(dl_lammps, "output_lammps.in")
In [ ]:
f = open("output_lammps.data")
print(f.read())
f.close()
In [ ]:
lammps_energy = dl_lammps.evaluate(utype="kcal * mol ** -1")
In [ ]:
# Compare energies
for k in energy_system1:
energy_difference = lammps_energy[k] - energy_system2[k]
print(k," difference:\t %.3f" % energy_difference)
Here, we demonstrate translating a small solvated peptide structure (from http://ambermd.org/tutorials/basic/tutorial0/index.htm) to LAMMPS using EEX.
In [ ]:
dl_dna = eex.datalayer.DataLayer("DNA_amber")
DNA_file = os.path.join("..", "examples", "amber","peptides", "alanine_dipeptide.prmtop")
eex.translators.amber.read_amber_file(dl_dna, DNA_file)
dl_dna.summary()
In [ ]:
eex.translators.lammps.write_lammps_file(dl_dna,"lammps_ala.data", "lammps_ala.in")
In [ ]:
f = open("lammps_ala.data")
print(f.read())
f.close()
In [ ]: